*----------------------------------------------------------------------------------------------------------	* 
* Silencing the Rails: A Study of the Noise-Safety Trade-off in Railway Quiet Zones                         *
* RESEARCHERS:		Emtiaz Hritan																		    *
* PROGRAMMED BY:	Emtiaz Hritan 					 											            *
* CREATED:			Oct. 25, 2023																		   	*
* LAST MODIFIED:	Oct. 7, 2025														       				*
*----------------------------------------------------------------------------------------------------------	*

clear all
set more off
* Set local paths
* Set this local datapath equal to the folder location for data
	local datapath "C:\Users\name\Replication Files Silencing the Rails\Data"
* Set this local outputpath equal to the folder location for outcome like tables
	local outputpath "C:\Users\name\Replication Files Silencing the Rails\Outcome"

* Install necessary packages
ssc install reghdfe
ssc install ftools
ssc install grstyle, replace
ssc install palettes, replace
ssc install colrspace, replace
ssc install ppmlhdfe, replace 
ssc install jwdid, replace 
ssc install csdid, replace 
ssc install drdid, replace 
ssc install did_imputation, replace
ssc install frause, replace
frause mpdta.dta, clear
ssc install hdfe, replace 
ssc install event_plot, replace
ssc install addplot, replace


*--------------------------------------------------------------------------
* Figure 4: Establishment of Quiet Zones in The U.S.
*--------------------------------------------------------------------------	
cd "`datapath'"
	
* Import the Uber/Lyft launch data
clear all
use quietzone
gen year=year(WhistleDate)
gen month =month(WhistleDate)
* Generate month of year variable
gen date_monthly = ym(year,month)
format %tm date_monthly
drop if WhistlebanCode==0
gen quietzone=1
collapse (sum) quietzone, by(year)

drop if year >2020
****Bar Graphs***
grstyle init
grstyle set plain
twoway (bar quietzone year, fcolor(none) lcolor(black)), ///
    xtitle("Year") ytitle("Number of Quiet Zones")
		
graph export "outputpath\Latex Source Files\quietzone_year.pdf", replace
*--------------------------------------------------------------------------
* Figure A5. :Establishment of Quiet Zones in The U.S. After 2005
*--------------------------------------------------------------------------	
* After 2005
	
drop if year ==2005
grstyle init
grstyle set plain
twoway (bar quietzone year, fcolor(none) lcolor(black)), ///
    xtitle("Year") ytitle("Number of Quiet Zones") ylabel(0(100)400)	
	

	
graph export "outputpath\Latex Source Files\quietzone_year_a2005.pdf", replace	
	
*--------------------------------------------------------------------------
* Figure 2a: Accidents at Railway Crossings From 1990 to 2019
*--------------------------------------------------------------------------	
* Let's see the overall accident trends 
* Seasonally adjusted
clear 
use final 
drop if year<1990
drop if year >2019
collapse (sum) accident, by(date_monthly year)
tsset date_monthly
tssmooth shwinters accident_sa = accident, additive from(0.3 0.1 0.1)


graph twoway (line accident_sa date_monthly, lcolor(blue)), ///
    ytitle("Number of Accidents") ///
    xtitle("Month") ///
    tline(2005m6, lcolor(green)) ///
    tline(1994m6, lcolor(green)) ///
    title("Accidents at Railway Crossings")
	
graph export "outputpath\Latex Source Files\trend_monthly_sa.pdf", replace	
*--------------------------------------------------------------------------
* Figure 2b: Accidents at Railway Crossings From 1990 to 2019
*--------------------------------------------------------------------------	
clear all 
use final 
drop if year<1990
drop if year >2019
drop if missing(accident)
collapse (sum) accident, by(year)


* Get 2010 base values by isolating 2005 and saving the sums in scalar variables
summarize accident if year == 2004
scalar base_accident = r(sum)



* Index each transportation mode by dividing by the 2005 base 
gen index_accident = (accident / base_accident) * 100


* Individual Indexed Trends
graph twoway ///
    (line index_accident year, lcolor(blue) lpattern(solid)) ///
    , title("") ///
    xtitle("Year") ytitle("Number of Accidents (Index, 2005=100)") xlabel(1990(5)2020)	
	
	
graph export "outputpath\Latex Source Files\trend_monthly_index.pdf", replace	
	
*--------------------------------------------------------------------------*
*  Figure A8. Panel (a): Yearly accident trends in top 10 cities by quiet zones intensity
*--------------------------------------------------------------------------*
***Create Quiet Zone Summary by City***

clear all
use final

* Clean
drop if year < 1990 | year > 2019
drop if missing(accident) | missing(cityname)

gen city_clean = lower(cityname)
gen city_std = cityname

replace city_std = "Chicago" if strpos(city_clean, "chicago") > 0
replace city_std = "New York City" if strpos(city_clean, "new york") > 0
replace city_std = "Los Angeles" if strpos(city_clean, "los angeles") > 0
replace city_std = "Houston" if strpos(city_clean, "houston") > 0
replace city_std = "Phoenix" if strpos(city_clean, "phoenix") > 0
replace city_std = "Philadelphia" if strpos(city_clean, "philadelphia") > 0
replace city_std = "San Antonio" if strpos(city_clean, "san antonio") > 0
replace city_std = "San Diego" if strpos(city_clean, "san diego") > 0
replace city_std = "Dallas" if strpos(city_clean, "dallas") > 0
replace city_std = "San Jose" if strpos(city_clean, "san jose") > 0

gen top10 = ///
  city_std == "New York City" | ///
  city_std == "Los Angeles"   | ///
  city_std == "Chicago"       | ///
  city_std == "Houston"       | ///
  city_std == "Phoenix"       | ///
  city_std == "Philadelphia"  | ///
  city_std == "San Antonio"   | ///
  city_std == "San Diego"     | ///
  city_std == "Dallas"        | ///
  city_std == "San Jose"


keep if top10 == 1
gen ever_qz = quietzone == 1
egen crossing_num = group(crossing_id)
* Create total crossings and quiet zone count per city
preserve
collapse (count) total_crossings=crossing_num (sum) num_qz=ever_qz, by(city_std)
gen qz_pct = 100 * num_qz / total_crossings
tempfile qz_stats
save `qz_stats'
restore


***Collapse Accident Counts and Merge Quiet Zone Stats***

collapse (sum) accident, by(city_std year)
merge m:1 city_std using `qz_stats', nogen


***Modify City Labels with Quiet Zone Stats***
   
*Create label for legend
* Use the following if want the intensity
*gen legend_label = city_std + " (" + string(round(qz_pct), "%2.0f") + "% QZ)"

gen legend_label = city_std  
  
* Sort to ensure line 1–10 corresponds to cities
sort city_std year

* Create macro labels (1 per city)
levelsof city_std, local(cities)
local i = 1
foreach c of local cities {
    preserve
    keep if city_std == "`c'"
    local lab`i' = legend_label[1]
    restore
    local ++i
}

* Now graph

twoway ///
  (line accident year if city_std == "Chicago", lcolor(black) lpattern(solid)) ///
  (line accident year if city_std == "New York City", lcolor(red) lpattern(shortdash)) ///
  (line accident year if city_std == "Los Angeles", lcolor(blue) lpattern(dot)) ///
  (line accident year if city_std == "Houston", lcolor(green) lpattern(longdash)) ///
  (line accident year if city_std == "Phoenix", lcolor(orange) lpattern(dash_dot)) ///
  (line accident year if city_std == "Philadelphia", lcolor(brown)) ///
  (line accident year if city_std == "San Antonio", lcolor(gray)) ///
  (line accident year if city_std == "San Diego", lcolor(purple)) ///
  (line accident year if city_std == "Dallas", lcolor(emerald)) ///
  (line accident year if city_std == "San Jose", lcolor(navy)) ///
  , ///
  legend(position(bottom) rows(3) ///
         label(1 "Chicago") ///
         label(2 "New York City") ///
         label(3 "Los Angeles") ///
         label(4 "Houston") ///
         label(5 "Phoenix") ///
         label(6 "Philadelphia") ///
         label(7 "San Antonio") ///
         label(8 "San Diego") ///
         label(9 "Dallas") ///
         label(10 "San Jose")) ///
  title("Accident Trends in Top 10 U.S. Cities") ///
  ytitle("Accidents") ///
  xtitle("Year") ///
  xline(2005, lpattern(solid) lcolor(black))  
  	
graph export "outputpath\Latex Source Files\top_10_city.pdf", replace

*--------------------------------------------------------------------------
*  Figure A8. Panel (b):  Yearly trends in accident counts across counties grouped by their quiet zone intensity
*--------------------------------------------------------------------------		

clear all 
use final 
drop if year<1990
drop if year >2019
drop if missing(accident)
	
gen ever_qz = quietzone == 1
gen county_fips= (statecode *1000) + countycode
egen crossing_num = group(crossing_id)

collapse (count) total_crossings=crossing_num (sum) num_qz=ever_qz, by(county_fips)	
gen qz_intensity = 100 * num_qz / total_crossings


tempfile quiet_zone_intensity
save `quiet_zone_intensity'	
	
* Get acident data 
clear 
use final 	
drop if year<1990
drop if year >2019
drop if missing(accident)
gen county_fips= (statecode *1000) + countycode
merge m:1 county_fips using `quiet_zone_intensity', force 	


xtile qz_group = qz_intensity if qz_intensity > 0, nq(4)
label define qzgrp 1 "Low" 2 "Medium" 3 "High" 4 "Very High"
label values qz_group qzgrp

collapse (sum) accident, by(qz_group year)

twoway (line accident year if qz_group==1, lcolor(blue) lpattern(dash)) ///
       (line accident year if qz_group==2, lcolor(black) lpattern(shortdash)) ///
       (line accident year if qz_group==3, lcolor(red) lpattern(solid)) ///
       (line accident year if qz_group==4, lcolor(orange) lpattern(longdash)) ///
       , legend(position(bottom) label(1 "Low QZ") label(2 "Med QZ") label(3 "High QZ") label(4 "Very High QZ")) ///
       title("Accident Trends by Quiet Zone Intensity Group") ///
       ytitle("Accidents (Count)") xtitle("Year") ///
  xline(2005, lpattern(solid) lcolor(black))
  
graph export "outputpath\Latex Source Files\quiet_zone_intensity.pdf", replace  
  
*--------------------------------------------------------------------------
*   Figure A4. : Accidents at Railway Crossings From 1975 to 2023
*--------------------------------------------------------------------------	  
  
* Let's see the overall accident trends 
clear all
use final 
collapse (sum) accident, by(date_monthly year)
* TWo-way graph with just count of accidents
graph twoway (line accident date_monthly, msymbol(I) msize(small) lcolor(blue)), ytitle(Number of Accidents) xtitle(Month) tline(2005m6, lcolor(green)) tline(1994m6, lcolor(green)) title(Accidents at Railway Crossings)  
  
  
graph export "jwdid_calender_never.eps", replace  
  
  
*--------------------------------------------------------------------------
*   Figure A6.: Kernel Density Plot of Monthly Accidents at Quiet Zone and Non-quiet Zone
*--------------------------------------------------------------------------
clear all
use final 

// Calculate the count of accidents per railway crossing per month
egen accidents_per_year = total(accident) if !missing(year), by(year)
egen accidents_per_month = total(accident) if !missing(date_monthly), by(date_monthly)

twoway (kdensity accidents_per_month if quietzone == 1, color(blue)) ///
       (kdensity accidents_per_month if quietzone == 0, color(red)), legend(label(1 "Quiet Zone") label(2 "Non-Quiet Zone")) ///
       title("Kernel Density Plot of Monthly Accidents") xtitle("Accidents per Month") ytitle("Density")  

graph export "outputpath\Latex Source Files\kernel_month.eps", replace	   
	   
*--------------------------------------------------------------------------
*   Figure A7. :Kernel Density Plot of Yearly Accidents at Quiet Zone and Non-quiet Zone
*--------------------------------------------------------------------------	   

// Create a kernel density plot for the entire dataset, using different colors/styles for each group

twoway (kdensity accidents_per_year if quietzone == 1, color(blue)) ///
       (kdensity accidents_per_year if quietzone == 0, color(red)), legend(label(1 "Quiet Zone") label(2 "Non-Quiet Zone")) ///
       title("Kernel Density Plot of Yearly Accidents") xtitle("Accidents per Year") ytitle("Density")
	     
  
graph export "outputpath\Latex Source Files\kernel_year.eps", replace  
  
  
********************************************************************************
* Figure A9. : Distribution of Train Volumes at Rail Crossings by Quiet Zone Status 
********************************************************************************

* Import and merge the all crossing, quiet zone and accident data 
clear all
use all_crossings
merge m:1 crossing_id using quietzone.dta, force 
replace quietzone=0 if quietzone==.
drop _merge

* Get one unique variable for the number of night and day train 
gen day_train_num=daythru 
replace day_train_num=TotalDaylightThruTrains if daythru==.

gen night_train_num=nghtthru
replace night_train_num=TotalNighttimeThruTrains if nghtthru==.

*How many trains travel on a track?

* Drop negative values
drop if night_train_num <0 
drop if day_train_num <0 

* Total train traffic
gen total_trains = day_train_num + night_train_num

*Trim observations below 1st and above 99th percentiles:
summarize total_trains, detail
gen dropflag = total_trains < r(p1) | total_trains > r(p99)
drop if dropflag
*Check distribution and summary by quiet zone status
tabstat day_train_num night_train_num total_trains, by(quietzone) stats(mean sd p50 min max n)



*****Panel (a)*****


*Kernel density by quiet zone status
drop if total_trains < 0 | total_trains > 250

twoway (kdensity total_trains if quietzone==1, lcolor(navy) lpattern(solid)) ///
       (kdensity total_trains if quietzone==0, lcolor(maroon) lpattern(dash)), ///
       legend(order(1 "Quiet Zone" 2 "Non-Quiet Zone") position(6)) ///
       title("Kernel Density of Train Traffic by Quiet Zone Status") ///
       xlabel(, grid) ylabel(, grid) ///
       xtitle("Total Trains per Day") ytitle("Density")

 graph export "outputpath\Latex Source Files\train_volume_1.pdf", as(pdf) name("Graph")

*****Panel (a)*****
	   
*Histogram by Quiet Zone Status
twoway ///
    (histogram total_trains if quietzone==1, percent color(navy%40) width(5) legend(label(1 "Quiet Zone"))) ///
    (histogram total_trains if quietzone==0, percent color(maroon%40) width(5) legend(label(2 "Non-Quiet Zone"))), ///
    title("Train Volume by Quiet Zone Status") ///
    legend(order(1 "Quiet Zone" 2 "Non-Quiet Zone") position(6)) ///
    ylabel(, grid) xlabel(0(10)80) ///
    xtitle("Total Trains per Day")


 graph export "outputpath\Latex Source Files\train_volume_4.pdf", as(pdf) name("Graph")	  
  
  
  
*--------------------------------------------------------------------------
* END
*--------------------------------------------------------------------------  
  
  
  
  
  
  